{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# elastic_constants_static calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The elastic_constants_static calculation style computes the elastic constants, $C_{ij}$, for a system by applying small strains and performing static energy minimizations of the initial and strained configurations.  Three estimates of the elastic constants are returned: one for applying positive strains, one for applying negative strains, and a normalized estimate that averages the &pm; strains and the symmetric components of the $C_{ij}$ tensor.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- Unlike the previous LAMMPS_ELASTIC calculation, this calculation does *not* perform a box relaxation on the system prior to evaluating the elastic constants.  This allows for the static elastic constants to be evaluated for systems that are relaxed to different pressures.\n",
    "- The elastic constants are estimated using small strains.  Depending on the potential, the values for the elastic constants may vary with the size of the strain.  This can come about either if the strain exceeds the linear elastic regime.\n",
    "- Some classical interatomic potentials have discontinuities in the fourth derivative of the energy function with respect to position.  If the strained states straddle one of these discontinuities the resulting static elastic constants values will be nonsense.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "The calculation method used here for computing elastic constants is based on the method used in the ELASTIC demonstration script created by Steve Plimpton and distributed with LAMMPS.\n",
    "\n",
    "The math in this section uses Voigt notation, where indicies i,j correspond to 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, and 6=xy, and x, y and z are orthogonal box vectors.\n",
    "\n",
    "A LAMMPS simulation performs thirteen energy/force minimizations\n",
    "\n",
    "- One for relaxing the initial system.\n",
    "\n",
    "- Twelve for relaxing systems in which a small strain of magnitude $\\Delta \\epsilon$ is applied to the system in both the positive and negative directions of the six Voigt strain components, $\\pm \\Delta \\epsilon_{i}$.\n",
    "\n",
    "The system virial pressures, $P_{i}$, are recorded for each of the thirteen relaxed configurations.  Two estimates for the $C_{ij}$ matrix for the system are obtained as\n",
    "\n",
    "$$ C_{ij}^+ = - \\frac{P_i(\\Delta \\epsilon_j) - P_i(0)}{\\Delta \\epsilon},$$\n",
    "\n",
    "$$ C_{ij}^- = - \\frac{P_i(0) - P_i(-\\Delta \\epsilon_j)}{\\Delta \\epsilon}.$$\n",
    "\n",
    "The negative out front comes from the fact that the system-wide stress state is $\\sigma_i = -P_i$.  A normalized, average estimate is also obtained by averaging the positive and negative strain estimates, as well as the symmetric components of the tensor\n",
    "\n",
    "$$ C_{ij} = \\frac{C_{ij}^+ + C_{ij}^- + C_{ji}^+ + C_{ji}^-}{4}.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2021-03-08 using iprPy version 0.10.4\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'elastic_constants_static'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is loaded from the database for the prototype."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc', database=library)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 108\n"
     ]
    }
   ],
   "source": [
    "sizemults = [3, 3, 3]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __strainrange__ specifies the $\\Delta \\epsilon$ strain range to use in estimating $C_{ij}$.\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "strainrange = 1e-7\n",
    "energytolerance = 1e-8\n",
    "forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. cij.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('cij.template', 'w') as f:\n",
    "    f.write(\"\"\"# Performs simulations to statically evaluate elastic constants using small strains\n",
    "# Based on the LAMMPS_ELASTIC script by Aidan Thompson (Sandia, athomps@sandia.gov)\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_info>\n",
    "\n",
    "change_box all triclinic\n",
    "\n",
    "# Specify strain\n",
    "variable strain equal <strainrange>\n",
    "\n",
    "# Define minimization parameters\n",
    "variable etol equal <etol>\n",
    "variable ftol equal <ftol>\n",
    "variable maxiter equal <maxiter>\n",
    "variable maxeval equal <maxeval>\n",
    "variable dmax equal <dmax>\n",
    "\n",
    "# Specify variables of the initial configuration's dimensions\n",
    "variable lx0 equal $(lx)\n",
    "variable ly0 equal $(ly)\n",
    "variable lz0 equal $(lz)\n",
    "\n",
    "# Specify the thermo properties to calculate\n",
    "variable peatom equal pe/atoms\n",
    "\n",
    "# Read in potential and thermo information\n",
    "include potential.in\n",
    "\n",
    "# Relax initial configuration and save as restart\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "write_restart initial.restart\n",
    "\n",
    "# Apply -xx strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${lx0}\n",
    "change_box all x delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +xx strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${lx0}\n",
    "change_box all x delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply -yy strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${ly0}\n",
    "change_box all y delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +yy strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${ly0}\n",
    "change_box all y delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply -zz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${lz0}\n",
    "change_box all z delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +zz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${lz0}\n",
    "change_box all z delta 0 ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply -yz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${lz0}\n",
    "change_box all yz delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +yz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${lz0}\n",
    "change_box all yz delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply -xz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${lz0}\n",
    "change_box all xz delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +xz strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${lz0}\n",
    "change_box all xz delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply -xy strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal -${strain}*${ly0}\n",
    "change_box all xy delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\n",
    "# Apply +xy strain\n",
    "clear\n",
    "read_restart initial.restart\n",
    "include potential.in\n",
    "\n",
    "variable delta equal ${strain}*${ly0}\n",
    "change_box all xy delta ${delta} remap units box\n",
    "minimize ${etol} ${ftol} ${maxiter} ${maxeval}\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. potential.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('potential.template', 'w') as f:\n",
    "    f.write(\"\"\"# NOTE: This script can be modified for different pair styles \n",
    "# See in.elastic for more info.\n",
    "\n",
    "# Choose potential\n",
    "<atomman_pair_info>\n",
    "\n",
    "# Setup minimization style\n",
    "min_modify dmax ${dmax}\n",
    "\n",
    "# Setup output\n",
    "thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe\n",
    "thermo_modify format float %.13e\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. elastic_constants_static()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def elastic_constants_static(lammps_command, system, potential, mpi_command=None,\n",
    "                             strainrange=1e-6, etol=0.0, ftol=0.0, maxiter=10000,\n",
    "                             maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Repeatedly runs the ELASTIC example distributed with LAMMPS until box\n",
    "    dimensions converge within a tolerance.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    strainrange : float, optional\n",
    "        The small strain value to apply when calculating the elastic\n",
    "        constants (default is 1e-6).\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is 10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'a_lat'** (*float*) - The relaxed a lattice constant.\n",
    "        - **'b_lat'** (*float*) - The relaxed b lattice constant.\n",
    "        - **'c_lat'** (*float*) - The relaxed c lattice constant.\n",
    "        - **'alpha_lat'** (*float*) - The alpha lattice angle.\n",
    "        - **'beta_lat'** (*float*) - The beta lattice angle.\n",
    "        - **'gamma_lat'** (*float*) - The gamma lattice angle.\n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the relaxed system.\n",
    "        - **'stress'** (*numpy.array*) - The measured stress state of the\n",
    "          relaxed system.\n",
    "        - **'C_elastic'** (*atomman.ElasticConstants*) - The relaxed system's\n",
    "          elastic constants.\n",
    "        - **'system_relaxed'** (*atomman.System*) - The relaxed system.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "\n",
    "    # Convert hexagonal cells to orthorhombic to avoid LAMMPS tilt issues\n",
    "    if am.tools.ishexagonal(system.box):\n",
    "        system = system.rotate([[2,-1,-1,0], [0, 1, -1, 0], [0,0,0,1]])\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    # Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='init.dat',\n",
    "                              units=potential.units,\n",
    "                              atom_style=potential.atom_style)\n",
    "    lammps_variables['atomman_system_info'] = system_info\n",
    "    lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)\n",
    "    lammps_variables['strainrange'] = strainrange\n",
    "    lammps_variables['etol'] = etol\n",
    "    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "    lammps_variables['maxiter'] = maxiter\n",
    "    lammps_variables['maxeval'] = maxeval\n",
    "    lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])\n",
    "    \n",
    "    # Fill in template files\n",
    "    template_file = 'cij.template'\n",
    "    lammps_script = 'cij.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "\n",
    "    template_file2 = 'potential.template'\n",
    "    lammps_script2 = 'potential.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file2, filedict)\n",
    "    with open(lammps_script2, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "    \n",
    "    # Run LAMMPS\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    \n",
    "    # Pull out initial state\n",
    "    thermo = output.simulations[0]['thermo']\n",
    "    pxx0 = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])\n",
    "    pyy0 = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])\n",
    "    pzz0 = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])\n",
    "    pyz0 = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])\n",
    "    pxz0 = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])\n",
    "    pxy0 = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])\n",
    "    \n",
    "    # Negative strains\n",
    "    cij_n = np.empty((6,6))\n",
    "    for i in range(6):\n",
    "        j = 1 + i * 2\n",
    "        # Pull out strained state\n",
    "        thermo = output.simulations[j]['thermo']\n",
    "        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])\n",
    "        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])\n",
    "        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])\n",
    "        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])\n",
    "        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])\n",
    "        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])\n",
    "        \n",
    "        # Calculate cij_n using stress changes\n",
    "        cij_n[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,\n",
    "                             pyz - pyz0, pxz - pxz0, pxy - pxy0]) / strainrange\n",
    "    \n",
    "    # Positive strains\n",
    "    cij_p = np.empty((6,6))\n",
    "    for i in range(6):\n",
    "        j = 2 + i * 2\n",
    "        # Pull out strained state\n",
    "        thermo = output.simulations[j]['thermo']\n",
    "        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])\n",
    "        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])\n",
    "        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])\n",
    "        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])\n",
    "        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])\n",
    "        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])\n",
    "        \n",
    "        # Calculate cij_p using stress changes\n",
    "        cij_p[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,\n",
    "                              pyz - pyz0, pxz - pxz0, pxy - pxy0]) / -strainrange\n",
    "    \n",
    "    # Average symmetric values\n",
    "    cij = (cij_n + cij_p) / 2\n",
    "    for i in range(6):\n",
    "        for j in range(i):\n",
    "            cij[i,j] = cij[j,i] = (cij[i,j] + cij[j,i]) / 2\n",
    "    \n",
    "    # Define results_dict\n",
    "    results_dict = {}\n",
    "    results_dict['raw_cij_negative'] = cij_n\n",
    "    results_dict['raw_cij_positive'] = cij_p\n",
    "    results_dict['C'] = am.ElasticConstants(Cij=cij)\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = elastic_constants_static(lammps_command, system, potential,\n",
    "                                        mpi_command = mpi_command,\n",
    "                                        strainrange = strainrange,\n",
    "                                        etol = energytolerance,\n",
    "                                        ftol = forcetolerance,\n",
    "                                        maxiter = maxiterations,\n",
    "                                        maxeval = maxevaluations,\n",
    "                                        dmax = maxatommotion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['raw_cij_negative', 'raw_cij_positive', 'C'])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __pressure_unit__ is the unit of pressure to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure_unit = 'GPa'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Print raw Cij values for $\\pm \\Delta \\epsilon$ states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raw Cij for negative strains (GPa) =\n",
      "[ 247.8622  147.8285  147.8285    0.0000    0.0000   -0.0000]\n",
      "[ 147.8285  247.8622  147.8285    0.0000    0.0000   -0.0000]\n",
      "[ 147.8285  147.8285  247.8622   -0.0000    0.0000   -0.0000]\n",
      "[  -0.0000    0.0001    0.0001  124.8381   -0.0000   -0.0000]\n",
      "[   0.0001   -0.0000    0.0001   -0.0000  124.8381   -0.0000]\n",
      "[   0.0001    0.0001   -0.0000   -0.0000   -0.0000  124.8381]\n",
      "\n",
      "Raw Cij for positive strains (GPa) =\n",
      "[ 247.8625  147.8283  147.8283   -0.0000    0.0000    0.0000]\n",
      "[ 147.8283  247.8625  147.8283    0.0000   -0.0000    0.0000]\n",
      "[ 147.8283  147.8283  247.8625    0.0000   -0.0000    0.0000]\n",
      "[   0.0000   -0.0001   -0.0001  124.8381    0.0000    0.0000]\n",
      "[  -0.0001    0.0000   -0.0001    0.0000  124.8381    0.0000]\n",
      "[  -0.0001   -0.0001    0.0000    0.0000    0.0000  124.8381]\n"
     ]
    }
   ],
   "source": [
    "print('Raw Cij for negative strains ('+pressure_unit+') =')\n",
    "for Ci in uc.get_in_units(results_dict['raw_cij_negative'], pressure_unit):\n",
    "    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))\n",
    "print()\n",
    "\n",
    "print('Raw Cij for positive strains ('+pressure_unit+') =')\n",
    "for Ci in uc.get_in_units(results_dict['raw_cij_positive'], pressure_unit):\n",
    "    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. Print averaged and refined Cij values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) =\n",
      "[ 247.8624  147.8284  147.8284    0.0000    0.0000    0.0000]\n",
      "[ 147.8284  247.8624  147.8284    0.0000    0.0000    0.0000]\n",
      "[ 147.8284  147.8284  247.8624    0.0000    0.0000    0.0000]\n",
      "[   0.0000    0.0000    0.0000  124.8381    0.0000    0.0000]\n",
      "[   0.0000    0.0000    0.0000    0.0000  124.8381    0.0000]\n",
      "[   0.0000    0.0000    0.0000    0.0000    0.0000  124.8381]\n"
     ]
    }
   ],
   "source": [
    "print('Cij ('+pressure_unit+') =')\n",
    "for Ci in uc.get_in_units(results_dict['C'].Cij, pressure_unit):\n",
    "    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
